Time delay is ubiquitous in most physical and biological systems, arising from finite propagation speeds of signals, which has been widely studied in the context of coupled limit cycle oscillator systems,^[1–4] for example, the semiconductor lasers,^[5] neuronal networks,^[6–8] etc.

Multiple myeloma is a hematologic malignancy that interferes in a destructive way with the natural bone renewal process,^[9–13] coupled with several diseases which involve abnormal and uncontrolled proliferation of tumor cells. Multiple myeloma-induced bone disease occurs due to dysregulated bone remodeling, where tumor cells interacting with the bone marrow microenvironment disrupt healthy kinetics, which is dependent on the behaviors of osteoclasts and osteoblasts, and ultimately leads to significant bone destruction.^[14–16] Currently, the therapy of multiple myeloma-induced bone disease is still a challenge by virtue of its incurability, but it may be controlled to some extent.^[17,18] Some mathematical models^[19–22] have been proposed to reveal the kinetics of bone cell populations in healthy and pathological bone remodeling, and underlying cellular and signaling mechanisms. It was demonstrated that the skeleton undergoes constant remodeling through bone resorption carried out by active osteoclasts and bone formation by osteoblasts, restoring the bone mass to its initial value, and exhibits regular periodic oscillations in bone cell populations. A limit cycle is observed in the phase plane of osteoclasts–osteoblasts.

Although there are a number of theoretical and experimental studies in various aspects of bone physiology, little work deals with the role of time delay in multiple myeloma-induced bone disease. However, there is increasing evidence for the existence of the delay during the bone remodeling cycles.^[23,24] For example, it is observed that the increase in bone cell populations takes place after 10–15 days in mice and 6–8 days in humans through therapy with G-CSF.^[24] Moreover, some clinical observations also demonstrated that signals transmission delays are inherent to bone systems due to the finite propagation speeds in the bone remodeling process.^[25–27]

Interesting questions now arise: how do the time delays influence the oscillation in a number of osteoclasts and osteoblasts? What are the effects of time delays on bone tumor (i.e., myeloma cells) in response to different treatments? To address these issues, according to the models^[19–22] of bone remodeling with multiple myeloma, we propose a mathematical model of bone remodeling kinetics with delay times of both osteoclast-derived paracrine signaling of tumor and tumor-derived paracrine signaling of osteoclast in this paper. Our model can characterize both the kinetics of bone remodeling during the disease progression and the time delays of osteoclasts on tumor cells due to the promotion of bone resorption for the growth of the tumor.

The rest of this paper is organized as follows. In Section 2, the mathematical model consisting of pathological bone remodeling and different therapeutic approaches is proposed. In Section 3, the influences of both osteoclast-derived paracrine signaling of tumor and tumor-derived paracrine signaling of osteoclast are discussed in multiple myeloma-induced bone disease, and the evolution of the bone system and the tumor in response to different treatments are studied in the presence of time delays. We draw some conclusions from the present study in Section 4.

In the presence of multiple myeloma, the bone remodeling process involves the activities of pre-osteoclasts, osteoclasts, pre-osteoblasts, osteoblasts, and myeloma cells. The scheme of interactions between bone cells and myeloma cells is shown by Fig. 1.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) Schematic diagram of interactions between bone cells and myeloma cells, showing (1) pre-osteoclasts differentiating into osteoclasts; (2) osteoclasts autocrine regulation; (3) osteoclast-derived factors stimulating osteoblasts production; (4) myeloma cell-derived factors affecting osteoclasts (PTHrP, TGF-β, IL-8, among others); (5) osteoclast-derived factors affecting myeloma cells; (6) pre-osteoblasts differentiating into osteoblasts; (7) osteoblasts autocrine regulation; (8) osteoblast-derived factors inhibiting osteoclasts production; (9) myeloma cells derived factors inhibiting osteoblasts production (DKK1, BMP, IGFs, among others).

Power law approximations are effective tools for describing and analyzing highly nonlinear biochemical systems.^[28] Bone remodeling is a nonlinear biological process involving complex mechanical and biochemical signaling pathways. Power law approximations are applied to the descriptions of biochemical autocrine and paracrine factors for bone system by Komarova et al., and are widely used in the subsequent studies of bone remodeling process.^[19–21] When the delay times of both osteoclast-derived paracrine signaling of tumor and tumor-derived paracrine signaling of osteoclast are considered, the time evolutions of osteoclasts, osteoblasts, and myeloma cells can be written as

(1)

(2)

(3)

(4)

On the other hand, for the multiple myeloma-induced bone disease, there are three possible pharmaceuticals for the treatment of disease: bisphosphonates which can suppress the production of osteoclasts,^[29–31] proteasome inhibitors which can stimulate osteoblast production and bone formation,^[16,32,33] and the anti-tumor therapy ^[13,20,29] which inhibits the growth of tumor cells. Thus, in the bone remodeling process involving the activities of osteoclasts, osteoblasts, and myeloma cells, the time evolution of cell populations Eqs. (1)–(3) can be adapted by

(5)

(6)

(7)

(8)

The above bone remodeling model can be studied by numerical simulation. Here, we take the parameters in accordance with Refs. [19] and [20]: , , , , , , , , , , , , , , , and , unless otherwise mentioned. The unit of time is days.

3.1. Oscillation phenomena of bone remodeling in the absence of time delays

In the absence of time delays, the kinetic behaviors of bone remodeling without or with myeloma cells are shown in Fig. 2. For normal bone remodeling, the system exhibits regular periodic oscillations in bone cell populations and a limit cycle is observed in the phase plane of osteoclasts–osteoblasts. In contrast, for pathological bone remodeling, the number of bone cells displays damped oscillations due to the profusion of myeloma cells, and the periods of oscillations become shorter than those in the normal case. It can be found that the amplitude of osteoclasts oscillations increases initially, reaches a maximum, and then decreases due to the growth of tumor with increasing time (see Figs. 2(a) and 2(c)). It means that the system tries to prevent bone resorption from abnormally coupling with bone formation, and tends to maintain normal bone remodeling. The oscillation amplitude of the number of osteoclasts increases after the invasion of myeloma cells, which is consistent with the experimental results.^[34,35] Figure 2(b) shows that the amplitude of osteoblasts oscillations remarkably decreases due to the invasion of tumor cells, which is a critical sign of multiple myeloma-induced bone disease. Figures 2(d)–2(f) show the variations in bone cell population; there is not a limit cycle in the phase plane of osteoclasts nor osteoblasts, and the oscillation amplitude of the ratio of bone cells C(t)/B(t) is reduced, leading to a reduction in bone mass. The variations in bone cell populations and bone mass of the mathematical model are in agreement with some experimental results,^[36–39] and similar results can be obtained in the following sections.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Kinetic behaviors of the system without time delays. The time evolution of the number of osteoclasts (a), the number of osteoblasts (b), the tumor (c), the ratio of C(t) to B(t) (e), and the bone mass (f). (d) The numbers of both osteoclasts and osteoblasts in phase plane.

3.2. Effects of time delay on osteoclast-derived paracrine signaling of tumor

3.2.1. The case of

The osteoclast-derived paracrine signaling of tumor can be characterized by the parameter . When is small, figures 3(a)–3(c) show that the growth of myeloma cells slows down, and the oscillation period of the ratio of osteoclasts to osteoblasts C(t)/B(t) is extended remarkably with the increase of time delay τ, leading to apparent difference in bone mass. When is large, however, the obvious alternations in the system by changing τ are not present, which indicates that the system is particularly sensitive to the variation of τ for small . In fact, although there is no time delay term in the equation (see Eq. (3)) of myeloma cells, the influence of the delay on the tumor is reflected in the curve of the tumor growth, implying that tumor cells are influenced by time delay effect of bone cells. Figures 3(c) and 3(f) show that the bone mass for small is higher than that for large at the same time (see the vertical dotted line): the smaller the value of , the less the bone loss in the bone system is, and the longer time delays can reduce bone loss and improve skeletal structure.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) The kinetics of the system with . (a)–(c) ; (d)–(e) .

3.2.2. The case of

When , figure 4 shows that the dependence of the systemic kinetics on time delay is similar to that in the case (Fig. 3), which indicates that the bone system exhibits strong resistance on account of the invasion of myeloma cells for long time delays when osteoclast-derived paracrine signaling of tumor is small (see the vertical dotted line).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Kinetic behaviors of the system with , with ((a)–(c)) and 0.5 ((d)–(f)).

When , however, it is found that the system does not display apparent variations by varying for small nor large as shown in Fig. 5. There is a competition between the delay and the osteoclast-derived paracrine signaling of tumor in the development of bone disease. In particular, for , the delay is predominant over the osteoclast-derived paracrine signaling of tumor for small , while the parameter is crucial for large . Accordingly, the growth of tumor cells and the progression of skeletal destruction slow down due to the effect of time delay for small whereas there are slight changes in the kinetic behaviors of the system for large . Moreover, it is observed that the larger the value of , the faster the tumor increases as it reaches a maximal value. In fact, with the growth of the tumor, more PTHrPs are secreted by tumor cells to promote the production of RANKL, leading to increased TGF-β by bone resorption of osteoclasts, further contributing to the tumor growth and resulting in a decline in bone mass. This suggests a vicious cycle imposed by tumor cells and osteoclasts,^[40–43] which induces bone destruction further and exacerbates bone disease.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Kinetic behaviors of the system with and ((a), (b)) and 0.5 ((c), (d)).

3.3. Effects of time delay on tumor-derived paracrine signaling of osteoclast

3.3.1. The case of

The tumor-derived paracrine signaling of osteoclast can be characterized by the parameter . Figure 6 shows typical time-dependent curves of the system for different values of . For large value of , the growth of tumor is slower and the oscillation period of C(t)/B(t) is longer than those in the case of small value of when . Therefore, increasing time delay can suppress the growth of tumor cells and the deterioration of bone destruction in the case of large value of (see the vertical dotted line).

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Kinetic behaviors of the system with , and ((a)–(c)) and 0.008 ((d), (e)).

3.3.2. The case of

When , it is found that from Fig. 7 the influence of time delay on the kinetics of the system is similar to that in the case (see Fig. 6). When , however, the time delay cannot induce obvious changes of the system for small value of nor value of large as shown in Fig. 8. The extents of skeletal destruction and the progression of bone disease can be changed by time delay for large value of (see the vertical dotted line), and the kinetic characteristics of the system depend strongly on the combined effects of the tumor-derived paracrine signaling and time delay. It is concluded that there is an optimal coupling between the tumor-derived paracrine signaling of osteoclasts and time delay by the above analysis during the progression of bone disease, which suppresses the tumor growth and the regression of bone disease.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Kinetic behaviors of the system with and ((a), (b)) and 0.008 ((c), (d)).

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Kinetic behaviors of the system with and ((a), (b)) and 0.008 ((c), (d)).

3.3.3. Therapeutic scenarios

The kinetic behaviors of myeloma cells and bone system in the presence of time delays are simulated for anti-tumor therapy, then together with either bisphosphonates or proteasome inhibitors, and numerical results are presented in Fig. 9. Different kinds of therapeutic strategies in our model start at time and stop at . For convenience, the time delays are set to be in the following.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Kinetic behaviors of the system with in response to combined treatment. Panels (a) and (b) show the anti-tumor and bisphosphonates ; panels (c) and (d) display the anti-tumor and proteasome inhibitors .

Figures 9(a) and 9(b) show how the combined therapy of anti-tumor and bisphosphonates affects kinetic properties of the system. The change of time delay can induce apparent variation of the system for drug . From Figs. 9(a) and 9(c), it is observed that a rapid reduction in the tumor growth takes places in the early stage of treatment. This is mainly due to the anti-tumor therapy , and almost independent of the influence of time delay during this time. However, the differences in the kinetics of the system are displayed by varying time delay in the later stage of tumor development. For drug , tumor cells exhibit relatively slow growth by increasing τ when treatment is stopped whereas the system does not exhibit remarkable variation for drug . Moreover, the cycles in the phase plane are changed with the variation of time delay τ for drug , which are represented by the curves in Fig. 9(b). This suggests that drug and long time delay τ can suppress tumor regrowth and bone disease relapse. By the above analysis of therapeutic scenarios the combined therapy of anti-tumor and bisphosphonates signifies the anti-tumor therapy in shrinkage of tumor but less effective in long-term suppression of tumor growth and prevention of tumor relapse. In fact, the bone system is capable of restoring itself to the healthy remodeling cycle after the effect of tumor cells has been completely eliminated. However, in either case, the tumor is usually never completely eliminated when the tumor reaches its limiting value, which may result in tumor regrowth after the treatment has been interrupted.

In this paper, we proposed a mathematical model of bone remodeling with time delays for multiple myeloma-induced bone disease. The effects of time delays on myeloma cells and bone system are analyzed from two representative systemic parameters: osteoclast-derived paracrine signaling of tumor and tumor-derived paracrine signaling of osteoclast . For osteoclast-derived paracrine signaling , the system exhibits diverse properties with time delays. It is shown that the time delays can apparently slow down the development of tumor cells and reduce the seriousness of the skeletal destruction for small value of when . Moreover, it is found that there is a vicious cycle between osteoclasts activity and tumor activity^[40–43] in such a way that the tumor development enhances bone resorption of osteoclasts, which in turn promotes the growth of tumor, leading to further bone destruction and severe bone disease. For tumor-derived paracrine signaling , the kinetic properties of the system are different from those for . In the cases of , the growth of tumor cells and the progression of bone disease are evidently slowed down by increasing time delay for large value of . It is found that there is an optimal coupling between the tumor-derived paracrine signaling of osteoclasts and time delay during the progression of bone diseases, which suppresses the tumor growth and the regression of bone disease.

Our model with time delays for multiple myeloma-induced bone disease not only can be used to analyze many qualitative properties of time delays and systemic parameters, but also might provide guidelines for probing the mechanism of bone destruction and exploring more effective therapeutic approaches for multiple myeloma-induced bone disease.